An exact Riemann Solver for multidimensional 
special relativist ic hydrodynamics 

Jose A. Pons 1 , Jose M— Marti 1 and Ewald Muller 2 

^epartament d'Astronomia i Astrofisica, U. Valencia, 46100 Burjassot, Spain 
qv | 2 MPI fur Astrophysik, Karl-Schwarzschild-Str. 1, 85748 Garching, Germany 

■*— > . e-mail: jose.a.pons@uv.es 

O: 

\Q • Abstract 

We have generalised the exact solution of the Riemann problem in special 

. ' relativistic hydrodynamics M for arbitrary tangential flow velocities. The 

►> . solution is obtained by solving the jump conditions across shocks plus an 

fvj ' ordinary differential equation arising from the self-similarity condition along 

\Q , rarefaction waves, in a similar way as in purely normal flow. The dependence 

of the solution on the tangential velocities is analysed. This solution has been 
^^ ' used to build up an exact Riemann solver implemented in a multidimensional 

i~l , relativistic (Godunov-type) hydro-code. 
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^ ! 1 Introduction 

The decay of a discontinuity separating two constant initial states (Riemann prob- 
;| ■ leni) has played a very important role in the development of numerical hydrody- 

"£i \ namic codes in classical (Newtonian) hydrodynamics after the pioneering work of 

Godunov |2j. Nowadays, most modern high-resolution shock-capturing methods 
^ | [pi are based on the exact or approximate solution of Riemann problems between 

adjacent numerical cells and the development of efficient Riemann solvers has be- 
come a research field in numerical analysis in its own (see, e.g., |14|| ). 
jrt ■ Riemann solvers began to be introduced in numerical relativistic hydrodynam- 

ics at the beginning of the nineties and, presently, the use of high-resolution shock- 
capturing methods based on Riemann solvers is considered as the best strategy to 
solve the equations of relativistic hydrodynamics in nuclear physics (heavy ion col- 
lisions) and astrophysics (stellar core collapse, supernova explosions, extragalactic 
jets, gamma-ray bursts). This fact has caused a rapid development of Riemann 
solvers for both special and general relativistic hydrodynamics p, 0] . 

The main idea behind the solution of a Riemann problem (defined by two 
constant initial states, L and R, at left and right of their common contact surface) 
is that the self-similarity of the flow through rarefaction waves and the Rankinc- 
Hugoniot relations across shocks allow one to connect the intermediate states /* 
(/ = L, R) with their corresponding initial states, /. The analytical solution of the 
Riemann problem in classical hydrodynamics [Q] rests on the fact that the normal 
velocity in the intermediate states, vf , can be written as a function of the pressure 



>< 



pi r in that state (and the flow conditions in state I). Thus, once pj r is known, i>™ 
and all other unknown state quantities of /* can be calculated. 

In the case of relativistic hydrodynamics the same procedure can be followed 
[pi p|, the major difference with classical hydrodynamics stemming from the role 
of tangential velocities. While in the classical case the decay of the initial dis- 
continuity does not depend on the tangential velocity (which is constant across 
shock waves and rarefactions), in relativistic calculations the components of the 
flow velocity are coupled through the presence of the Lorentz factor. 

2 The equations of relativistic hydrodynamics 

The equations of relativistic hydrodynamics admit a conservative formulation 
which has been exploited in the last decade to implement high-resolution shock- 
capturing methods. In Minkowski space time the equations in this formulation 
read 

9 t U + 9iPW=0 (1) 

where U and F^(U) (i = 1,2,3) are, respectively, the vectors of conserved vari- 
ables and fluxes 

\J=(D,S\S 2 ,S\t) t (2) 

PW = {Dv\ S 1 v i +p6 li ,S 2 v i +p8 2i ,S 3 v i +p8 M ,S i -Dv i ) T - (3) 

The conserved variables (the rest-mass density, D, the momentum density, S l , 
and the energy density r) are defined in terms of the primitive variables,(p, v l , e), 
according to 

D = pW , S i = phW 2 v l , r = phW 2 -p-D (4) 

where W = (1 — v 2 ) -1 / 2 is the Lorentz factor and h = 1 + e + p/p the specific 
enthalpy. 

In the following wc shall restrict our discussion to an ideal gas equation of state 
with constant adiabatic exponent, 7, for which the specific internal energy is given 

by 

(7 - 1)P 

3 Relation between the normal flow velocity and 
pressure behind relativistic rarefaction waves 

Choosing the surface of discontinuity to be normal to the x-axis, rarefaction waves 
are self-similar solutions of the flow equations depending only on the combination 
£ = x/t. Getting rid of all the terms with y and z derivatives in equations (|lj) and 
substituting the derivatives of x and t in terms of the derivatives of £, the system 
of equations can be reduced to just one ordinary differential equation (ODE) and 



two algebraic conditions 



phW 2 (v x - £)dv x + (1 - iv x )dp = 




(G) 


hWv y = constant 




(7) 


hWv z = constant, 




(8) 


by 






v x (l - cj) ± c s V(l - v 2 )[l - v 2 c 2 s - {v x ) 2 {l - 


-c*)} 


(9) 


1 - v 2 c 2 


J 



with £ constrained by 

because non-trivial similarity solutions exist only if the Wronskian of the original 
system vanish. We have denoted by c s the speed of sound, provided by the equation 
of state. The plus and minus sign correspond to rarefaction waves propagating to 
the right (72—*) and left (??.<_), respectively. The two solutions for £ correspond 
to the maximum and minimum eigenvalues of the Jacobian matrix associated to 
jpO)(u) |H, p|, generalizing the result found for vanishing tangential velocity Q. 

From equations (0) and (g) it follows that v y /v z = constant, i.e., the tangential 
velocity does not change direction across rarefaction waves. Notice that, in a 
kinematical sense, the Newtonian limit (v l <C 1) leads to W = 1, but equations (Q) 
and (jg) do not reduce to the classical limit v y ' z — constant, because the specific 
enthalpy still couples the tangential velocities. Thus, even for slow flows, the 
Riemann solution presented in this paper must be employed for thermodynamically 
relativistic situations (h 3> 1). The same result can be deduced from the Rankinc- 
Hugoniot relations for shock waves (see next section) . 

Using (||) and the definition of the sound speed, the ODE (J6|) can be written 
as 

^! = ±^ x (io) 

dp W 2 1P v'l + stefc, "*,«*) 



where u* = \/{v y ) 2 + (v z ) 2 is the absolute value of the tangential velocity and 

(i -£ ±v x y 

Considering that in a Riemann problem the state ahead of the rarefaction wave 
is known, equation (|10|) can be integrated with the constraint hWv 1 = constant, 
allowing to connect the states ahead (a) and behind (6) the rarefaction wave. The 
solution is only a function of p>\, and can be stated in compact form as 

v* b =K<L( Pb ). (12) 

4 Relation between post-shock flow velocities and 
pressure for relativistic shock waves. 

The Rankine-Hugoniot conditions relate the states on both sides of a shock and 
are based on the continuity of the mass flux and the energy-momentum flux across 
shocks. Their relativistic version was first obtained by Taub pT| (see also O. pj). 



Considering the surface of discontinuity as normal to the x-axis, the invariant 
mass flux across the shock can be written as 



j = W s D a {V s - v x ) = W s D b {V s - v x b ). 



(13) 



where V s is the coordinate velocity of the hyper-surface that defines the position 
of the shock wave and W s is the correspondent Lorentz factor, W s = (1 — V 2 ) 1 / 2 . 
According to our definition, j is positive for shocks propagating to the right. In 
terms of the mass flux, j, the Rankinc-Hugoniot conditions are 



M 



j 




"1" 


w s 


D 


' 


-gx- 


s 


[d\ 



J_ 
w t 

[hWv v ] = 0, 
[hWv z ] = 0, 



[v p\ 



J 

W s 



D 



(14) 



(15) 

(16) 
(17) 

(18) 



Equations ( [Tq ) and ( |lT|) imply that the quantity hWv v ' z is constant across a shock 
wave and, hence, that the orientation of the tangential velocity does not change. 
The latter result also holds for rarefaction waves (see §3). Equations (fl4|), ( |15| ) 
and ( |l8| ) can be manipulated to obtain vf as a function of p b , j and V s . Using the 
relation S x = (r + p + D)v x and after some algebra, one finds 



, W s {p b -p a )\ ( s(W 8 vl 

y b = I KW a v x + ^- h a W a + (p b - Pa ) < 



3 



J 



PaW a 



(19) 

The final step is to express j and V s as a function of the post-shock pressure. 
First, from the definition of the mass flux we obtain 



V d 



pjWlv x a ± \jWf + p\Wl{\ - vf) 
P 2 a W*+f 



(20) 



where V+ (V~) corresponds to shocks propagating to the right (left). 

Second, from the Rankine-Hugoniot relations and the physical solution of h b 
obtained from the Taub adiabat [pL3|] (the relativistic version of the Hugoniot adi- 
abat), that relates only thermodynamic quantities on both sides of the shock, the 
square of the mass flux j can be obtained as a function of p b . Using the pos- 
itive (negative) root of j 2 for shock waves propagating towards the right (left), 
the desired relation between the post-shock normal velocity v b and the post-shock 
pressure p b is obtained. In a compact way the relation reads 

v x b =S^(p b ). (21) 

We refer to the interested reader to references M, ffl for further details. 



5 The solution of the Riemann problem with ar- 
bitrary tangential velocities. 

The time evolution of a Riemann problem can be represented as: 

/ -» LW^L*CR*W^R (22) 

where W denotes a simple wave (shock or rarefaction), moving towards the initial 
left (<— ) or right (— >) states. Between them, two new states appear, namely L* 
and i?», separated from each other through the third wave C, which is a contact 
discontinuity. Across the contact discontinuity pressure and normal velocity are 
constant, while the density and the tangential velocity exhibits a jump. 

As in the Newtonian case, the compressive character of shock waves (density 
and pressure rise across the shock) allows us to discriminate between shocks (<S) 
and rarefaction waves (TZ): 

WL M = { ^ M ' Pb ~ Pa (23) 

1 ' 1 <$<-(-►) . Pb>Pa 

where p is the pressure and subscripts a and b denote quantities ahead and behind 
the wave. For the Riemann problem a = L(R) and b = L*{R*) for WV_ and W_^, 
respectively. 

The solution of the Riemann problem consists in finding the positions of the 
waves separating the four states and the intermediate states, L* and R*. The 
functions W_+ and WV_ allow one to determine the functions v#*(p) and vf jlt (p), 
respectively. The pressure p* and the flow velocity v^ in the intermediate states 
are then given by the condition 



TV* 



(P*)=»L(P.)=<- (24) 



which is an implicit algebraic equation in p* and can be solved by means of an 
iterative method. When p* and v% have been obtained, the equation of state gives 
the specific internal energy and the remaining state variables of the intermediate 
state /* can be calculated using the relations between /* and the respective initial 
state / given through the corresponding wave. 

Notice that the solution of the Riemann problem depends on the modulus of 
v*, but not on the direction of the tangential velocity. Figure |] shows the solution 
of a particular Riemann problem jlQ] for different values of the tangential velocity 
v y = 0,0.5,0.9,0.99. The crossing point of any two lines in the upper panel gives 
the pressure and the normal velocity in the intermediate states. Whereas the 
pressure in the intermediate state can take any value between pl and pr, the 
normal flow velocity can be arbitrarily close to zero in the case of an extremely 
relativistic tangential flow. 

To study the influence of tangential velocities on the solution a Riemann prob- 
lem, we have calculated the solution of a standard test involving the propagation of 
a relativistic blast wave produced by a large jump in the initial pressure distribu- 
tion for different combinations of tangential velocities H . Although the structure 
of the solution remains unchanged for different tangential velocities, the values in 
the constant state may change by a large amount. 
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Figure 1: Graphical solution in the (p, v x ) plane (upper panel) and in the (y , v x )- 
plane (lower panel) of the relativistic Riemann problem with initial data pl = 1.0, 
pL = 1.0, u£ = 0.0; p_R = 0.1, pn = 0.125 and v R = 0.0 for different values of 
the tangential velocity v l = 0, 0.5, 0.9, 0.999, represented by solid, dashed, dashed- 
dotted and dotted lines, respectively. An ideal gas EOS with 7 = 1.4 was assumed. 
The crossing point of any two lines in the upper panel gives the pressure and the 
normal velocity in the intermediate states. The value of the tangential velocity in 
the states L* and R* is obtained from the value of the corresponding functions at 
v x in the lower panel, and Iq gives the solution for vanishing tangential velocity. 
The range of possible solutions is given by the shaded region in the upper panel. 




Figure 2: Exact (solid lines) and numerical profiles along the diagonal of pressure 
(squares) , density (crosses) and normal velocity (diamonds) for the 2-dimensional 
relativistic shock tube discussed in the text. 

6 Discussion and conclusions. 



We have obtained the exact solution of the Riemann problem in special relativistic 
hydrodynamics with arbitrary tangential velocities. Unlike in Newtonian hydro- 
dynamics, tangential velocities are coupled with the rest of variables through the 
Lorentz factor, present in all terms in all equations. It strongly affects the so- 
lution, especially for ultra-relativistic tangential flows. In addition, the specific 
enthalpy also acts as a coupling factor and modifies the solution for the tangential 
velocities in thermodynamically relativistic situations (energy density and pressure 
comparable to or larger than the proper rest-mass density) , rendering the classical 
solution incorrect in slow flows with very large internal energies. 

Our solution has interesting practical applications. First, it can be used to test 
the different approximate relativistic Riemann solvers and the multi-dimensional 
hydrodynamic codes based on directional splitting. Second, it can be used to con- 
struct multi-dimensional relativistic Godunov-type hydro codes. As an example, 
we have simulated a relativistic tube pQ| , in a 100 x 100 Cartesian grid, where 
the initial discontinuity was located in a main diagonal. The initial states were 
Pl = 10, pn, = 1, pl = 13.3, pr — 0.66E — 3, vl = 0, vr — 0, and the adiabatic 
index is 7 = 5/3. Spatial order of accuracy was set to second order by means of 
a monotonic piecewise linear reconstruction procedure and second order in time 
is obtained by using a Runge-Kutta method for time advancing. The exact solu- 
tion of the Riemann problem is used at every interface to calculate the numerical 
fluxes. The results are shown in Figure 2, and are comparable to those obtained 
with other HRSC methods. Profiles of all variables are stable and discontinuities 



are well resolved without excessive smearing. An efficient implementation of this 
exact Riemann solver in the context of multidimensional relativistic PPM m is in 
progress and will be reported elsewhere. 
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